## this code runs the analyses presented in the main text and supplemental information, 
## for all respondents who passed both attention checks. It should be run third. 

#### clear workspace,  load R libraries ####
rm(list=ls())
library(car)
library(ggplot2)
library(reshape2)
library(plyr)
library(survey)
library(tidyverse)
library(arm)
library(Cairo)
library(ggthemes)
library(modelr)
library(sparsereg)
library(viridis)
library(ggpubr)
library(stargazer)
library(broom)
library(xtable)
library(grid)
library(gridExtra)
library(RColorBrewer)
library(lmtest)

## set working directory here
setwd('~/Dropbox (Yale_FES)/ypccc_parrish/GND/replication')
setwd('~/Dropbox/GND')

vcovCluster <- function(
  model,
  cluster
)
{
  require(sandwich)
  require(lmtest)
  if(nrow(model.matrix(model))!=length(cluster)){
    stop("check your data: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)           
  K <- model$rank   
  if(M<50){
    warning("Fewer than 50 clusters, variances may be unreliable (could try block bootstrap instead).")
  }
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}

#### Attribute effects ####
#### read in pre-saved data ####
stack1<-read.csv('gnd_stack1.csv')
# levels(stack1$Economic)
stack1<-stack1%>%
  dplyr::mutate(Economic=relevel(Economic,ref="none"),
                Social=relevel(Social,ref="none"),
                Sponsor=relevel(Sponsor,ref="Democrats"))
stack2<-read.csv('gnd_stack2.csv')

#### function for 1st conjoint analysis ####
baselines<-data.frame('term'=c("Economic","Social","Carbon","Cost","Size","Sponsor",
                               "no economic policy","no social policy","no carbon tax","$10 per month","$100 billion per year",
                               "Democrats"),
                      'estimate'=0,'std.error'=0,'statistic'=0,'p.value'=0,'conf.low'=0,
                      'conf.high'=0,'model'="Conjoint 1",'type'=c(rep("label",6),rep("baseline",6)))
baselines$att<-rep(baselines$term[1:6],2)

## function to use weighted OLS 
lmfunc1<-function(x,y,z){## x is data; y is name of group column; z is value in group column
  l<-lm(selected~Economic+Social+Carbon+Cost+Size+Sponsor+factor(vignette_id),
        data=x,weights=weight)
  c<-coeftest(l,vcov = vcovCluster(l,factor(x[,"id"])))
  pevec <- c[,1] 
  sevec <- c[,2]
  out<-tidy(l,conf.int=TRUE)
  out$estimate<-pevec
  out$std.error<-sevec
  # out<-tidy(l,conf.int=TRUE)
  out<-out%>%
    dplyr::mutate(conf.low=estimate-(1.96*std.error),
                  conf.high=estimate+(1.96*std.error),
                  term=gsub("Economic","Economic: ",term),
                  term=gsub("Social","Social: ",term),
                  term=gsub("Carbon","Carbon: ",term),
                  term=gsub("Cost","Cost: ",term),
                  term=gsub("Size","Size: ",term),
                  term=gsub("Sponsor","Sponsor: ",term))%>%
    dplyr::mutate(att=gsub(":.*","",term))%>%
    dplyr::mutate(term=gsub("Economic: ","",term),
                  term=gsub("Social: ","",term),
                  term=gsub("Carbon: ","",term),
                  term=gsub("Cost: ","",term),
                  term=gsub("Size: ","",term),
                  term=gsub("Sponsor: ","",term))
  
  out$model<-"Conjoint 1"
  out$type<-"estimate"
  out<-rbind(baselines,out)
  out$term<-factor(out$term,levels=c("(Intercept)",
                                     "factor(vignette_id)2",
                                     "factor(vignette_id)3",
                                     "bipartisan",
                                     "Democrats",
                                     "Sponsor",
                                     "$55 per month",
                                     "$35 per month",
                                     "$10 per month",
                                     "Cost",
                                     "$500 billion per year",
                                     "$250 billion per year",
                                     "$100 billion per year",
                                     "Size",
                                     "revenue neutral tax",
                                     "tax and dividend",
                                     "tax and invest",
                                     "no carbon tax",
                                     "Carbon",
                                     "unionized clean energy jobs",
                                     "retrain fossil fuel workers",
                                     "job guarantee",
                                     "no economic policy",
                                     "Economic",
                                     "free college",
                                     "$15 minimum wage",
                                     "health insurance",
                                     "affordable housing",
                                     "no social policy",
                                     "Social"))
  out[,y]<-z  
  return(out)
  outsg<<-out%>%dplyr::filter(type!="label"&type!="baseline")%>%
    dplyr::select(term,estimate,std.error,conf.low,conf.high)
  
}


out1<-lmfunc1(stack1,"party","All")

sg<-function(mod,groupname){
  out<-mod%>%dplyr::filter(type!="label"&type!="baseline")
  out<-out[,c("term","estimate","std.error","conf.low","conf.high",groupname)]
  names(out)<-c("Term","Estimate","Std. Error","95 CI (H)","95 CI (L)",groupname)
  return(out)
}
out1sg<-sg(out1,"party")

stargazer(out1sg%>%dplyr::select(-c(party)),summary=FALSE,rownames=FALSE)

## extract t-values: 
summary(lm(selected~Economic+Social+Carbon+Cost+Size+Sponsor+factor(vignette_id),
data=stack1,weights=weight))


#### function for 2nd conjoint analysis ####
baselines2<-data.frame('term'=c("Transportation","Fossil Infrastructure","Carbon","Cost","Size","Investments","Energy",
                               "no transportation policy","no fossil infrastructure policy","no carbon tax","$10 per month","$100 billion per year",
                               "no investments","no energy policy"),
                      'estimate'=0,'std.error'=0,'statistic'=0,'p.value'=0,'conf.low'=0,
                      'conf.high'=0,'model'="Conjoint 1",'type'=c(rep("label",7),rep("baseline",7)))
baselines2$att<-rep(baselines2$term[1:7],2)
stack2<-stack2%>%
  dplyr::mutate(Legacy=relevel(Legacy,ref="no fossil infrastructure policy"),
                Energy=relevel(Energy,ref="no energy policy"),
                Investments=relevel(Investments,ref="no investments"),
                Transportation=relevel(Transportation,ref="no transportation policy"))

#### function for 2nd conjoint #### 
lmfunc2<-function(x,y,z){
  lm2<-lm(selected~Carbon+Legacy+Energy+Investments+Cost+Size+Transportation+factor(vignette_id),
          data=x,weights=weight)
  c<-coeftest(lm2,vcov = vcovCluster(lm2,factor(x[,"id"])))
  pevec <- c[,1] 
  sevec <- c[,2]
  out<-tidy(lm2,conf.int=TRUE)%>%dplyr::mutate(model="Conjoint 2")
  out$estimate<-pevec
  out$std.error<-sevec
  out<-out%>%
    dplyr::mutate(conf.low=estimate-(1.96*std.error),
                  conf.high=estimate+(1.96*std.error),
                  term=gsub("Transportation","Transportation: ",term),
                  term=gsub("Legacy","Fossil Infrastructure: ",term),
                  term=gsub("Carbon","Carbon: ",term),
                  term=gsub("Cost","Cost: ",term),
                  term=gsub("Size","Size: ",term),
                  term=gsub("Investments","Investments: ",term),
                  term=gsub("Energy","Energy: ",term))%>%
    dplyr::mutate(att=gsub(":.*","",term))
  out[out$term=="Investments: Carbon:  removal","term"]<-"Investments: direct air capture"
  out[out$term=="Carbon: tax and invest: ",'term']<-"Carbon: tax and invest"
  out<-out%>%
    dplyr::mutate(term=gsub("Transportation: ","",term),
                  term=gsub("Fossil Infrastructure: ","",term),
                  term=gsub("Carbon: ","",term),
                  term=gsub("Cost: ","",term),
                  term=gsub("Size: ","",term),
                  term=gsub("Investments: ","",term),
                  term=gsub("Energy: ","",term))
  
  out$term
  out$model<-"Conjoint 2"
  out$type<-"estimate"
  out<-rbind(baselines2,out)
  out$term<-factor(out$term,levels=c("(Intercept)",
                                       "factor(vignette_id)2",
                                       "factor(vignette_id)3",
                                     "$55 per month",
                                     "$35 per month",
                                     "$10 per month",
                                     "Cost",
                                     "$500 billion per year",
                                     "$250 billion per year",
                                     "$100 billion per year",
                                     "Size",
                                       "public transit",
                                       "electric car mandate",
                                       "no transportation policy",
                                       "Transportation",
                                       "prosecute",
                                       "eliminate subsidies",
                                     "coal plant securitization",
                                       "no fossil infrastructure policy",
                                       "Fossil Infrastructure",
                                       "natural carbon storage",
                                       "direct air capture",
                                       "building retrofits",
                                       "no investments",
                                       "Investments",
                                       "mandate coal phaseout 2030",
                                       "CES with nuclear",
                                       "CES with CCS",
                                       "CES renewables only",
                                       "no energy policy",
                                       "Energy",
                                       "revenue neutral tax",
                                       "tax and dividend",
                                       "tax and invest",
                                       "no carbon tax",
                                       "Carbon")) 
  out[,y]<-z
  return(out)
}
out2<-lmfunc2(stack2,"party","All")
out2sg<-sg(out2,"party")
stargazer(out2sg%>%dplyr::select(-c(party)),summary=FALSE,rownames=FALSE)

g_legend<-function(a.gplot){
  tmp <- ggplot_gtable(ggplot_build(a.gplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)}


#### main effects plot: conjoint 1 ####
m<-ggplot(out1[out1$term!="(Intercept)"&
                   grepl("vignette",out1$term)==FALSE,],aes(x=term,y=estimate,shape=type,alpha=type,color=att))+
  geom_pointrange(aes(ymin=conf.low,ymax=conf.high),
                  position=position_dodge(0.75),size=.4)+
  scale_shape_manual(values=c(1,19,20),guide=FALSE)+
  scale_alpha_manual(values=c("label"=0,"baseline"=1,"estimate"=1),guide=FALSE)+
  scale_x_discrete(position="top")+
  geom_hline(yintercept=0,lty="dashed",color="darkgray")+
  scale_y_continuous(limits=c(-.3,0.3),breaks=c(-.3,0,0.3))+
  labs(x=element_blank(),y=element_blank())+
  coord_flip()+
  theme_bw()+
  theme(axis.text=element_text(size=12),
        legend.text=element_text(size=12),
        axis.ticks.y=element_blank(),
        axis.text.y=element_blank(),
        legend.title=element_blank(),
        aspect.ratio=3,
        legend.position="none",
        plot.margin = unit(c(0, 0, 0, 0), "cm")
)
m


#### main effects plot: conjoint 2 ####
m2<-ggplot(out2[out2$term!="(Intercept)"&
                  grepl("vignette",out2$term)==FALSE,],aes(x=term,y=estimate,shape=type,alpha=type,color=att))+
  geom_pointrange(aes(ymin=conf.low,ymax=conf.high),
                  position=position_dodge(0.75),size=.4)+
  scale_shape_manual(values=c(1,19,20),guide=FALSE)+
  scale_alpha_manual(values=c("label"=0,"baseline"=1,"estimate"=1),guide=FALSE)+
  geom_hline(yintercept=0,lty="dashed",color="darkgray")+
  scale_y_continuous(limits=c(-.3,0.3),breaks=c(-.3,0,0.3))+
  scale_x_discrete(position="top")+
  labs(x=element_blank(),y=element_blank())+
  coord_flip()+
  theme_bw()+
  theme(axis.text=element_text(size=12),
        legend.text=element_text(size=12),
        axis.ticks.y=element_blank(),
        axis.text.y=element_blank(),
        legend.title=element_blank(),
        aspect.ratio=3,
        legend.position="none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"))
m2


#### split sample by party ####
rep<-stack1%>%dplyr::filter(PID=="Republican")
dem<-stack1%>%dplyr::filter(PID=="Democrat")
outr<-lmfunc1(rep,"party","Republicans")
outd<-lmfunc1(dem,"party","Democrats")
out.party<-do.call(rbind,list(outd,outr))
out.party$party<-factor(out.party$party,
                        levels=c("Democrats","Republicans","All"))
pcolors<-c("Democrats"="blue","Republicans"="red","All"="darkgray")
out.party$party[out.party$type=="baseline"]="All"
stargazer(sg(out.party,"party"),summary=FALSE,rownames=FALSE)
## extract t stats 
out.party

#### fig function for conjoint 1 ####
c1plot<-function(x,y,z,t){ ## x is data; y is grouping var for color;
  # z is vector of colors; t is breaks for color scale
  y=sym(y)
  x<-x[x[,'term']!="(Intercept)"&
         grepl("vignette",x[,'term'])==FALSE,]
  m<-ggplot(x,aes(x=term,y=estimate,alpha=type,shape=type,color=!!y))+
    geom_pointrange(aes(ymin=conf.low,ymax=conf.high),
                    position=position_dodge(.75),size=.4)+
    scale_shape_manual(values=c(1,19,20),guide=FALSE)+
    scale_alpha_manual(values=c("label"=0,"baseline"=1,"estimate"=1),guide=FALSE)+
    scale_color_manual(values=z, breaks=t)+
    scale_x_discrete(position="bottom")+
    geom_hline(yintercept=0,lty="dashed",color="darkgray")+
    scale_y_continuous(limits=c(-.35,0.35),breaks=c(-.3,0,0.3))+
    labs(x=element_blank(),y=element_blank())+
    coord_flip()+
    theme_bw()+
    theme(axis.text=element_text(size=12),
          legend.text=element_text(size=12),
          axis.ticks.y=element_blank(),
          axis.text.y=element_text(hjust=0.5,face=c(rep("plain",2),"bold",
                                                    rep("plain",3),"bold",
                                                    rep("plain",3),"bold",
                                                    rep("plain",4),"bold",
                                                    rep("plain",4),"bold",
                                                    rep("plain",5),"bold")),
          legend.title=element_blank(),
          aspect.ratio=3,
          legend.position="bottom",
          plot.margin = unit(c(0, 0, 0, -3.8), "cm"))
  return(m)
}

#### C1 plot without labels ####
c1plot_nolabels<-function(x,y,z,t){ ## x is data; y is grouping var for color;
  # z is vector of colors; t is breaks for color scale
  y=sym(y)
  x<-x[x[,'term']!="(Intercept)"&
         grepl("vignette",x[,'term'])==FALSE,]
  m<-ggplot(x,aes(x=term,y=estimate,alpha=type,shape=type,color=!!y))+
    geom_pointrange(aes(ymin=conf.low,ymax=conf.high),
                    position=position_dodge(0.75),size=.4)+
    scale_shape_manual(values=c(1,19,20),guide=FALSE)+
    scale_alpha_manual(values=c("label"=0,"baseline"=1,"estimate"=1),guide=FALSE)+
    scale_color_manual(values=z, breaks=t)+
    scale_x_discrete(position="bottom")+
    geom_hline(yintercept=0,lty="dashed",color="darkgray")+
    scale_y_continuous(limits=c(-.35,0.35),breaks=c(-.3,0,0.3))+
    labs(x=element_blank(),y=element_blank())+
    coord_flip()+
    theme_bw()+
    theme(axis.text=element_text(size=12),
          legend.text=element_text(size=12),
          axis.ticks.y=element_blank(),
          axis.text.y=element_blank(),
          legend.title=element_blank(),
          aspect.ratio=3,
          legend.position="bottom",
          plot.margin = unit(c(0, 0, 0, -3.8), "cm"))
  return(m)
}
dr<-c1plot(out.party,"party",pcolors,c("Democrats", "Republicans"))
dr

party.leg <- g_legend(dr)

dr <- dr + theme(legend.position="none")

quartz(width=7, height=5.5)
grid.arrange(m, dr, nrow=1)

writefig<-function(name,...){
  png(file=paste0("figures/",name,".png"),width=7,height=5.5,units="in",#width=720,height=1080,units="px",
      # width=6,height=4,units='in',
      res=400)
}
writefig("conjoint1-netpartisan")
ggarrange(m, dr, nrow=1)
dev.off()




#### conjoint 2: split sample by party ####
dem2<-stack2%>%dplyr::filter(PID=="Democrat")
rep2<-stack2%>%dplyr::filter(PID=="Republican")
outd2<-lmfunc2(dem2,"party","Democrats")
outr2<-lmfunc2(rep2,"party","Republicans")
out.party2<-do.call(rbind,list(outd2,outr2))
out.party2$party<-factor(out.party2$party,
                        levels=c("Democrats","Republicans","All"))

pcolors<-c("Democrats"="blue","Republicans"="red","All"="darkgray")
out.party2$party[out.party2$type=="baseline"]<-"All"
stargazer(sg(out.party2,"party"),summary=FALSE,rownames=FALSE)

#### fig function for conjoint 2 ####
c2plot<-function(x,y,z,t){
  y=sym(y)
  x<-x[x[,'term']!="(Intercept)"&
         grepl("vignette",x[,'term'])==FALSE,]
  m<-ggplot(x,aes(x=term,y=estimate,alpha=type,shape=type,color=!!y))+
    geom_pointrange(aes(ymin=conf.low,ymax=conf.high),
                    position=position_dodge(0.75),size=.4)+
    scale_shape_manual(values=c(1,19,20),guide=FALSE)+
    scale_alpha_manual(values=c("label"=0,"baseline"=1,"estimate"=1),guide=FALSE)+
    scale_color_manual(values=z,breaks=t)+
    scale_x_discrete(position="bottom")+
    geom_hline(yintercept=0,lty="dashed",color="darkgray")+
    scale_y_continuous(limits=c(-.35,0.35),breaks=c(-.3,0,0.3))+
    labs(x=element_blank(),y=element_blank())+
    coord_flip()+
    theme_bw()+
    theme(axis.text=element_text(size=12),
          legend.text=element_text(size=12),
          axis.ticks.y=element_blank(),
          axis.text.y=element_text(hjust=0.5,face=c(rep("plain",3),"bold",
                                                    rep("plain",3),"bold",
                                                    rep("plain",3),"bold",
                                                    rep("plain",4),"bold",
                                                    rep("plain",4),"bold",
                                                    rep("plain",5),"bold",
                                                    rep("plain",4),"bold")),
          legend.title=element_blank(),
          aspect.ratio=3,
          legend.position="bottom",
          plot.margin = unit(c(0, 0, 0, -3.8), "cm"))
  return(m)
}
c2plot_nolabels<-function(x,y,z,t){
  y=sym(y)
  x<-x[x[,'term']!="(Intercept)"&
         grepl("vignette",x[,'term'])==FALSE,]
  m<-ggplot(x,aes(x=term,y=estimate,alpha=type,shape=type,color=!!y))+
    geom_pointrange(aes(ymin=conf.low,ymax=conf.high),
                    position=position_dodge(0.75),size=.4)+
    scale_shape_manual(values=c(1,19,20),guide=FALSE)+
    scale_alpha_manual(values=c("label"=0,"baseline"=1,"estimate"=1),guide=FALSE)+
    scale_color_manual(values=z,breaks=t)+
    scale_x_discrete(position="bottom")+
    geom_hline(yintercept=0,lty="dashed",color="darkgray")+
    scale_y_continuous(limits=c(-.35,0.35),breaks=c(-.3,0,0.3))+
    labs(x=element_blank(),y=element_blank())+
    coord_flip()+
    theme_bw()+
    theme(axis.text=element_text(size=12),
          legend.text=element_text(size=12),
          axis.ticks.y=element_blank(),
          axis.text.y=element_blank(),
          legend.title=element_blank(),
          aspect.ratio=3,
          legend.position="bottom",
          plot.margin = unit(c(0, 0, 0,0), "cm"))
  return(m)
}
dr2<-c2plot(out.party2,"party",pcolors,c("Democrats","Republicans"))
dr2
dr2 <- dr2 + theme(legend.position="none")

quartz(width=7, height=5.5)
grid.arrange(m2, dr2, nrow=1)
writefig("conjoint2-netpartisan")
ggarrange(m2, dr2, nrow=1)
dev.off()



#### conjoint 1: split sample by race ####
white<-stack1%>%dplyr::filter(race_Q=="White/Caucasian")
black<-stack1%>%dplyr::filter(race_Q=="Black/African American")
hisp<-stack1%>%dplyr::filter(race_Q=="Hispanic/Latino")
other<-stack1%>%dplyr::filter(race_Q=="Asian"|race_Q=="Other"|race_Q=="Native American/Pacific Islander")


lmw<-lmfunc1(white,"race","White")
lmb<-lmfunc1(black,"race","Black")
lmh<-lmfunc1(hisp,"race","Hispanic")
lmo<-lmfunc1(other,"race","Other")
names(out1)[11]<-"race"
out.race<-do.call(rbind,list(lmw,lmb,lmh))
out.race$race[out.race$type=="baseline"]<-"All"

rcolors<-c(brewer.pal(9,"Set1")[9],brewer.pal(12,"Paired")[c(1,8,10)])

stargazer(sg(out.race,"race"),rownames=FALSE,summary=FALSE)


### use this function to create conjoint 1 figures without y axis labels
r1plot<-function(x,y,z,t){ ## x is data; y is grouping var for color; z is vector of colors; 
  # t is vector of breaks for color scale
  y=sym(y)
  x<-x[x[,'term']!="(Intercept)"&
         grepl("vignette",x[,'term'])==FALSE,]
  m<-ggplot(x,aes(x=term,y=estimate,alpha=type,shape=type,color=!!y))+
    geom_pointrange(aes(ymin=conf.low,ymax=conf.high),
                    position=position_dodge(.75),size=.4)+
    scale_shape_manual(values=c(1,19,20),guide=FALSE)+
    scale_alpha_manual(values=c("label"=0,"baseline"=1,"estimate"=1),guide=FALSE)+
    scale_color_manual(values=z, breaks=t)+
    scale_x_discrete(position="bottom")+
    geom_hline(yintercept=0,lty="dashed",color="darkgray")+
    scale_y_continuous(limits=c(-.45,0.45),breaks=c(-.3,0,0.3))+
    labs(x=element_blank(),y=element_blank())+
    coord_flip()+
    theme_bw()+
    theme(axis.text=element_text(size=12),
          legend.text=element_text(size=12),
          axis.ticks.y=element_blank(),
          axis.text.y=element_blank(),
          legend.title=element_blank(),
          aspect.ratio=3,
          legend.position="bottom",
          plot.margin = unit(c(0, 0, 0, 0), "cm"))
  return(m)
}
drace<-r1plot(out.race,"race",rcolors,c("White", "Black","Hispanic"))
drace




#### conjoint 1: split sample by income ####
low<-stack1%>%dplyr::filter(income=="<$50K")
med<-stack1%>%dplyr::filter(income=="$50K-$100K")
high<-stack1%>%dplyr::filter(income=="$100K+")
lmlow<-lmfunc1(low,"income","<$50K")
lmmed<-lmfunc1(med,"income","$50K-$100K")
lmhigh<-lmfunc1(high,"income","$100K+")
out.income<-out1
names(out.income)[11]<-"income"

out.income<-do.call(rbind,list(lmlow,lmmed,lmhigh))#,out.income))
out.income$income<-factor(out.income$income,levels=c("All",
                                                     "<$50K","$50K-$100K","$100K+"))
out.income$income[out.income$type=="baseline"]<-"All"
icolors<-rcolors
stargazer(sg(out.income,"income"),rownames=FALSE,summary=FALSE)

i<-c1plot(out.income,"income",icolors,c("<$50K","$50K-$100K","$100K+"))
i

quartz(width=7, height=5.5)
grid.arrange(drace, i, nrow=1)
writefig("conjoint1-subgroup")
ggarrange(drace,i,nrow=1)
dev.off()





#### conjoint 1: split sample by rural/urban ####
rural <-stack1%>%dplyr::filter(metro=="rural")
urban<-stack1%>%dplyr::filter(metro=="urban")

lmrural<-lmfunc1(rural,"metro","rural")
lmurban<-lmfunc1(urban,"metro","urban")
out.metro<-out1
names(out.metro)[11]<-"metro"

out.metro<-do.call(rbind,list(lmrural, lmurban))#,out.income))
out.metro$metro<-factor(out.metro$metro,levels=c("All",
                                                     "rural","urban"))
out.metro$metro[out.metro$type=="baseline"]<-"All"
mcolors<-rcolors
stargazer(sg(out.metro,"metro"),rownames=FALSE,summary=FALSE)

m<-c1plot(out.metro,"metro",mcolors,c("rural","urban"))
m

quartz(width=7, height=5.5)
grid.arrange(m, nrow=1)
dev.copy2pdf(file="figures/conjoint1-urban-rural.pdf")





## conjoint2
rural2<-stack2%>%dplyr::filter(metro=="rural")
urban2<-stack2%>%dplyr::filter(metro=="urban")

lmrural2<-lmfunc2(rural2,"metro","rural")
lmurban2<-lmfunc2(urban2,"metro","urban")
out.metro2<-out2
names(out.metro2)[11]<-"metro"

out.metro2<-do.call(rbind,list(lmrural2, lmurban2))
out.metro2$metro[out.metro2$type=="baseline"]<-"All"

stargazer(sg(out.metro2,"metro"),rownames=FALSE,summary=FALSE)


m2<-c2plot(out.metro2,"metro",mcolors,c("rural", "urban"))
m2

quartz(width=7, height=5.5)
grid.arrange(m2, nrow=1)
dev.copy2pdf(file="figures/conjoint2-urban-rural.pdf")

#### conjoint 2: race ####
#### conjoint 1: split sample by race ####
white<-stack2%>%dplyr::filter(race_Q=="White/Caucasian")
black<-stack2%>%dplyr::filter(race_Q=="Black/African American")
hisp<-stack2%>%dplyr::filter(race_Q=="Hispanic/Latino")
other<-stack2%>%dplyr::filter(race_Q=="Asian"|race_Q=="Other"|race_Q=="Native American/Pacific Islander")


lmw<-lmfunc2(white,"race","White")
lmb<-lmfunc2(black,"race","Black")
lmh<-lmfunc2(hisp,"race","Hispanic")
lmo<-lmfunc2(other,"race","Other")
names(out2)[11]<-"race"
out.race2<-do.call(rbind,list(lmw,lmb,lmh))
out.race2$race[out.race2$type=="baseline"]<-"All"
stargazer(sg(out.race2,"race"),rownames=FALSE,summary=FALSE)
drace2<-c2plot_nolabels(out.race2,"race",rcolors,c("White", "Black","Hispanic"))


#### conjoint 2: income ####
low2<-stack2%>%dplyr::filter(income=="<$50K")
med2<-stack2%>%dplyr::filter(income=="$50K-$100K")
high2<-stack2%>%dplyr::filter(income=="$100K+")
lmlow2<-lmfunc2(low2,"income","<$50K")
lmmed2<-lmfunc2(med2,"income","$50K-$100K")
lmhigh2<-lmfunc2(high2,"income","$100K+")
out.income2<-out2
names(out.income2)[11]<-"income"
out.income2<-do.call(rbind,list(lmlow2,lmmed2,lmhigh2))
out.income2$income<-factor(out.income2$income,levels=c("All",
                                                     "<$50K","$50K-$100K","$100K+"))
out.income2$income[out.income2$type=="baseline"]<-"All"
stargazer(sg(out.income2,"income")[1:44,],rownames=FALSE,summary=FALSE)
stargazer(sg(out.income2,"income")[45:nrow(sg(out.income2,"income")),],rownames=FALSE,summary=FALSE)


i2 <- c2plot(out.income2,"income",icolors,c("<$50K","$50K-$100K","$100K+"))


quartz(width=7, height=5.5)
grid.arrange(drace2, i2, nrow=1)
writefig("conjoint2-subgroup")
ggarrange(drace2,i2,nrow=1)
dev.off()

#### Conjoint 1: Additional AMCEs ####
# Age
table(stack1$age_Q)
age1<-stack1%>%dplyr::filter(age_Q=="18-24")
age2<-stack1%>%dplyr::filter(age_Q=="25-34")
age3<-stack1%>%dplyr::filter(age_Q=="35-44")
age4<-stack1%>%dplyr::filter(age_Q=="45-54")
age5<-stack1%>%dplyr::filter(age_Q=="55-64")
age6<-stack1%>%dplyr::filter(age_Q=="65+")
lma1<-lmfunc1(age1,"age","18-24")
lma2<-lmfunc1(age2,"age","25-34")
lma3<-lmfunc1(age3,"age","35-44")
lma4<-lmfunc1(age4,"age","45-54")
lma5<-lmfunc1(age5,"age","55-64")
lma6<-lmfunc1(age6,"age","65+")
out.age<-do.call(rbind,list(lma1,lma2,lma3,lma4,lma5,lma6))
out.age$age<-factor(out.age$age,levels=c("All","18-24","25-34",
                                         "35-44","45-54","55-64","65+"))
out.age$age[out.age$type=="baseline"]<-"All"
acolors<-c(brewer.pal(9,"Set1")[9],brewer.pal(12,"Paired")[c(1,3,5,8,10,12)])

nrow(sg(out.age,"age"))
stargazer(sg(out.age,"age")[1:44,],rownames=FALSE,summary=FALSE)
stargazer(sg(out.age,"age")[45:88,],rownames=FALSE,summary=FALSE)
stargazer(sg(out.age,"age")[89:nrow(sg(out.age,"age")),],rownames=FALSE,summary=FALSE)


a<-c1plot(out.age,"age",acolors,
          c("18-24","25-34",
            "35-44","45-54","55-64","65+"))
quartz(width=7, height=5.5)
a
writefig("supp_c1_age")
a
dev.off()

# Gender
table(stack1$gender_Q)
g1<-stack1%>%dplyr::filter(gender_Q=="Female")
g2<-stack1%>%dplyr::filter(gender_Q=="Male")
lmg1<-lmfunc1(g1,"gender","Female")
lmg2<-lmfunc1(g2,"gender","Male")
out.gender<-do.call(rbind,list(lmg1,lmg2))
out.gender$gender<-factor(out.gender$gender,levels=c("All","Female","Male"))
out.gender$gender[out.gender$type=="baseline"]<-"All"
gcolors<-c(brewer.pal(8,"Set2")[8],brewer.pal(12,"Paired")[c(8,10)])
# gcolors<-c("All"="darkgray","Female"=rainbow_hcl(2)[1],"Male"=rainbow_hcl(2)[2])
stargazer(sg(out.gender,"gender"),rownames=FALSE,summary=FALSE)
# 
# 
g<-c1plot(out.gender,"gender",gcolors,
          c("Female","Male"))
quartz(width=7, height=5.5)
g
writefig("supp_c1_gender")
g
dev.off()

# Census Region
table(stack1$region_name)
r1<-stack1%>%dplyr::filter(region_name=="Midwest Region")
r2<-stack1%>%dplyr::filter(region_name=="Northeast Region")
r3<-stack1%>%dplyr::filter(region_name=="South Region")
r4<-stack1%>%dplyr::filter(region_name=="West Region")
lmr1<-lmfunc1(r1,"region","Midwest")
lmr2<-lmfunc1(r2,"region","Northeast")
lmr3<-lmfunc1(r3,"region","South")
lmr4<-lmfunc1(r4,"region","West")
out.region<-do.call(rbind,list(lmr1,lmr2,lmr3,lmr4))
out.region$region<-factor(out.region$region,levels=c("All","Midwest","Northeast","South","West"))
out.region$region[out.region$type=="baseline"]<-"All"
regcolors<-c(brewer.pal(8,"Set2")[8],brewer.pal(12,"Paired")[c(1,5,8,10)])
stargazer(sg(out.region,"region")[1:44,],rownames=FALSE,summary=FALSE)
stargazer(sg(out.region,"region")[45:nrow(sg(out.region,"region")),],rownames=FALSE,summary=FALSE)

# 
# 
r<-c1plot(out.region,"region",regcolors,
          c("Midwest","Northeast","South","West"))
quartz(width=7, height=5.5)
r
writefig("supp_c1_region")
r
dev.off()

## conjoint 2
# Age
age1<-stack2%>%dplyr::filter(age_Q=="18-24")
age2<-stack2%>%dplyr::filter(age_Q=="25-34")
age3<-stack2%>%dplyr::filter(age_Q=="35-44")
age4<-stack2%>%dplyr::filter(age_Q=="45-54")
age5<-stack2%>%dplyr::filter(age_Q=="55-64")
age6<-stack2%>%dplyr::filter(age_Q=="65+")
lma1<-lmfunc2(age1,"age","18-24")
lma2<-lmfunc2(age2,"age","25-34")
lma3<-lmfunc2(age3,"age","35-44")
lma4<-lmfunc2(age4,"age","45-54")
lma5<-lmfunc2(age5,"age","55-64")
lma6<-lmfunc2(age6,"age","65+")
out.age<-do.call(rbind,list(lma1,lma2,lma3,lma4,lma5,lma6))
out.age$age<-factor(out.age$age,levels=c("All","18-24","25-34",
                                         "35-44","45-54","55-64","65+"))
out.age$age[out.age$type=="baseline"]<-"All"
# 
nrow(sg(out.age,"age"))
stargazer(sg(out.age,"age")[1:66,],rownames=FALSE,summary=FALSE)
stargazer(sg(out.age,"age")[67:nrow(sg(out.age,"age")),],rownames=FALSE,summary=FALSE)


a<-c2plot(out.age,"age",acolors,
          c("18-24","25-34",
            "35-44","45-54","55-64","65+"))
quartz(width=7, height=5.5)
a
writefig("supp_c2_age")
a
dev.off()
# 
# # Gender
g1<-stack2%>%dplyr::filter(gender_Q=="Female")
g2<-stack2%>%dplyr::filter(gender_Q=="Male")
lmg1<-lmfunc2(g1,"gender","Female")
lmg2<-lmfunc2(g2,"gender","Male")
out.gender<-do.call(rbind,list(lmg1,lmg2))
out.gender$gender<-factor(out.gender$gender,levels=c("All","Female","Male"))
out.gender$gender[out.gender$type=="baseline"]<-"All"
stargazer(sg(out.gender,"gender"),rownames=FALSE,summary=FALSE)
#
#
g<-c2plot(out.gender,"gender",gcolors,
          c("Female","Male"))
quartz(width=7, height=5.5)
g
writefig("supp_c2_gender")
g
dev.off()
# 
# # Census Region
r1<-stack2%>%dplyr::filter(region_name=="Midwest Region")
r2<-stack2%>%dplyr::filter(region_name=="Northeast Region")
r3<-stack2%>%dplyr::filter(region_name=="South Region")
r4<-stack2%>%dplyr::filter(region_name=="West Region")
lmr1<-lmfunc2(r1,"region","Midwest")
lmr2<-lmfunc2(r2,"region","Northeast")
lmr3<-lmfunc2(r3,"region","South")
lmr4<-lmfunc2(r4,"region","West")
out.region<-do.call(rbind,list(lmr1,lmr2,lmr3,lmr4))
out.region$region<-factor(out.region$region,levels=c("All","Midwest","Northeast","South","West"))
out.region$region[out.region$type=="baseline"]<-"All"
stargazer(sg(out.region,"region"),rownames=FALSE,summary=FALSE)
#
#
r<-c2plot(out.region,"region",regcolors,
          c("Midwest","Northeast","South","West"))
quartz(width=7, height=5.5)
r
writefig("supp_c2_region")
r
dev.off()


## median response time
below<-stack1%>%dplyr::filter(Duration..in.seconds.>=828)
above<-stack1%>%dplyr::filter(Duration..in.seconds.<828)
lmabove<-lmfunc1(above,"responsetime","Above")
lmbelow<-lmfunc1(below,"responsetime","Below")
out.med<-do.call(rbind,list(lmabove,lmbelow))
out.med$responsetime<-factor(out.med$responsetime,levels=c("All","Above","Below"))
out.med$responsetime[out.med$type=="baseline"]<-"All"
mcolors<-gcolors
stargazer(sg(out.med,"responsetime"),rownames=FALSE,summary=FALSE)
# 
# 
m<-c1plot(out.med,"responsetime",mcolors,
          c("Above","Below"))
quartz(width=7, height=5.5)
m
writefig("supp_c1_responsetime")
m
dev.off()

## median response time
below2<-stack2%>%dplyr::filter(Duration..in.seconds.>=828)
above2<-stack2%>%dplyr::filter(Duration..in.seconds.<828)
lmabove2<-lmfunc2(above2,"responsetime","Above")
lmbelow2<-lmfunc2(below2,"responsetime","Below")
out.med2<-do.call(rbind,list(lmabove2,lmbelow2))
out.med2$responsetime<-factor(out.med2$responsetime,levels=c("All","Above","Below"))
out.med2$responsetime[out.med2$type=="baseline"]<-"All"
stargazer(sg(out.med2,"responsetime"),rownames=FALSE,summary=FALSE)
# 
# 
m2<-c2plot(out.med2,"responsetime",mcolors,
          c("Above","Below"))
quartz(width=7, height=5.5)
m2
writefig("supp_c2_responsetime")
m2
dev.off()

